unit StochasticProcess;
interface
uses ODEModel, ODEIntegrator, Math, Helper, SysUtils;

type

TWhiteNoise = class
private
  variance : extended;
  values : Vector;
public
  constructor Create(passBand : extended);
  function Get(t : extended ) : extended;
end;

TStochasticProcess = class(TODEModel)
public
  BandPass : extended;
  CorrelationWindow : extended;
  function Get(n : integer) : extended;  virtual; abstract;
protected
  Tau, _T, K, Etta : extended;
  noise : TWhiteNoise;
end;

// Process1
TStochasticProcess1 = class(TStochasticProcess)
public
  constructor Create(integrator: TODEIntegrator; D, L : extended);
  procedure F(t: extended; y:Vector; var dx : Vector); override;
  function Get(n : integer) : extended; override;
end;

// Process2
TStochasticProcess2 = class(TStochasticProcess)
public
  constructor Create(integrator: TODEIntegrator; D, L, w0 : extended);
  procedure F(t: extended; y:Vector; var dx : Vector); override;
  function Get(n : integer) : extended; override;
end;

implementation

  constructor TWhiteNoise.Create(passBand : extended);
  begin
    variance := passBand/(Pi*2);
  end;

  function TWhiteNoise.Get(t : extended) : extended;
  var
    position : integer;
  begin

    position := Trunc(t/(1/variance));

    if( position > High(values) ) then
    begin
      SetLength(values, Length(values) +1);
      values[High(values)] := Randg(0.0, sqrt(variance));
    end;

    result := values[position];
  end;

{ stochastic 2}
  constructor TStochasticProcess2.Create(integrator: TODEIntegrator; D, L, w0 : extended);
  var b0, b2, a4, a2, a0 : extended;
  begin
    problemSize := 2;
    SetLength(InitialConditions,2);

    InitialConditions[0] := 0;
    InitialConditions[1] := 0;
    CorrelationWindow := L/(sqr(L) + sqr(w0));

    K := sqrt(2*D*L/(sqr(L) + sqr(w0)));
    Tau := sqrt(1/(sqr(L)+sqr(w0)));
    _T := Tau;
    Etta := L/sqrt(sqr(L) + sqr(w0));

    BandPass := 100/(_T);

    Inherited Create(integrator);
  end;

  procedure TStochasticProcess2.F(t: extended; y:Vector; var dx : Vector);
  var
    x : extended;
  begin
    if( noise = nil ) then
        noise := TWhiteNoise.Create(BandPass);
      x := noise.Get(t);
      dx[0] := (K*x-y[1])/sqr(_T) - 2*Etta*y[0]/_T;
      dx[1] := y[0];
  end;

  function TStochasticProcess2.Get(n : integer) : extended;
  begin
    result := Tau*y[n, 0] + y[n, 1];
  end;

{ stochastic 1 }
  constructor TStochasticProcess1.Create(integrator: TODEIntegrator; D, L : extended);
  begin
    problemSize := 2;
    SetLength(InitialConditions,2);
    InitialConditions[0] := 0;
    InitialConditions[1] := 0;
    CorrelationWindow := 1/L;
    self.K := sqrt(2*D/L);
    _T := 1/L;
    BandPass := 100/_T;


    Inherited Create(integrator);
  end;

  function TStochasticProcess1.Get(n : integer) : extended;
  begin
    result := Y[n, 0];
  end;

  procedure TStochasticProcess1.F(t: extended; y:Vector; var dx : Vector);
    var
      x: extended;
    begin
      if( noise = nil ) then
        noise := TWhiteNoise.Create(BandPass);
      x := noise.Get(t);
      dx[0] := (K*x - y[0])/_T;
    end;
end.
